########################################################
###### Replication Code for Figure 4 and Table A7 ######
################## January 9th, 2025 ################### 
########################################################

rm(list=ls())

### Script produces Table A1 and saves it in `/tables` folder

# load libraries
library(ggplot2)
library(dplyr)
library(modelsummary)
library(kableExtra)
library(stringr)
library(fixest)
library(stargazer)
library(gridExtra)
library(flextable)

# load main data
df <- readRDS("data/workfile.rds")

########### Descriptives
desc_table <- df %>% 
  ungroup() %>%
  filter(!duplicated(meta_UUID)) %>%
  select(dem_female, dem_age, dem_rural, dem_university, lrscale,
         mean_outparty_feeling, d_affiliated, high_income,
         medium_income, low_income, income_missing, 
         trust_in_opposing_parties, political_interest_num) %>% 
  mutate_all(as.numeric)

##### Declare covariates ####
vars <- colnames(desc_table)[colnames(desc_table) != "meta_Weight"] # remove weights



##### Declare covariates ####
treats <- c("control", "treated_econ", "treated_culture")
vars <- c("dem_female", "dem_age", "dem_rural", "lrscale", "d_affiliated", 
          "high_income", "medium_income", "low_income", "income_missing",
          "dem_university", "mean_outparty_feeling", "mean_close_num", 
          "we_vs_they_party", "high_pol_interest")

##### Define helper function #####
# get_stats returns mean and sd for defined covariates (vars) and subgroups (treats)
get_stats <- function(var, treats, data) {
  c(mean(data[, var][data[, treats] ==1 ], na.rm = T), sd(data[, var][data[, treats] ==1 ], na.rm = T))
}


## loop
btable <- lapply(treats, function(y) t(sapply(vars, function(x) {
  get_stats(x, y, data = df)
})))

## to df, add n, colnames
btable <- data.frame(do.call(cbind, btable))

## times 100
perc <- c("dem_rural", "dem_female", "dem_university", "high_income", 
          "medium_income", "low_income", "income_missing", "high_pol_interest",
          "d_affiliated")
btable[rownames(btable) %in% perc, ] <- btable[rownames(btable) %in% perc, ] * 100

## add varnames
btable$variable <- rownames(btable)

## add n obs
n <- sapply(vars, function(x) sum(!is.na(df[, x])))
btable <- cbind(n, btable)


## change colnames
cn <- paste0(rep(treats, each = 2), rep(c("_mean", "_sd"), 2))
colnames(btable)[c(-1, -ncol(btable))] <- cn
btable[, -which(colnames(btable) == "variable")] <- 
  data.frame(round(btable[, -which(colnames(btable) == "variable")], 1))
btable$Variable <- btable$variable %>% 
  recode("dem_female" = "Female (in %)", 
         "dem_age" = "Age", 
         "dem_rural" = "Rural (in %)", 
         "dem_university" = "University Education (in %)", 
         "close_num" = "Social distance", 
         "feeling_toward" = "Feeling score", 
         "lrscale" = "Left-right scale", 
         "high_income" = "High income (in %)", 
         "medium_income" = "Medium income (in %)",
         "low_income" = "Low income (in %)",
         "income_missing" = "Missing income (in %)",
         "mean_outparty_feeling" = "Average Out-Party Feeling", 
         "mean_close_num" = "Average Out-Party Social Distance", 
          "we_vs_they_party" = "Use of we vs. they", 
         "high_pol_interest" = "High political interest (in %)",
         "d_affiliated" = "Partisans (in %)"
         )

rownames(btable) <- NULL
btable <- btable[, c(ncol(btable), 2:(ncol(btable)-2))]

##### Latex output via stargazer
stargazer(btable, summary =F, type = "html", out = "tables/TabA1_desc_table.html", rownames = F, digits = 1)
?stargazer
